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Abstract 



^ . A spectral integral method (lEM) for solving the two-body Schrodinger equation in configura- 

^ I tion space is generalized to the calculation of the corresponding T— matrix. It is found that the 

P^ ' desirable features of the lEM, such as the economy of mesh-points for a given required accuracy, 

are carried over also to the solution of the T— matrix. However the algorithm is considerably more 
^~-*', complex, because the T— matrix is a function of two variables r and r', rather than only one vari- 

able r, and has a slope discontinuity at r = r'. For a simple exponential potential an accuracy of 
7 significant figures is achieved, with the number N of Chebyshev support points in each partition 
y^ . equal to 17. For a potential with a large repulsive core, such as the potential between two He 

£^ ' atoms, the accuracy decreases to 4 significant figures, but is restored to 7 if A^ is increased to 65. 
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I. INTRODUCTION 

The two-body T— matrix is most often used for the calculation of two-body scattering 
phase shifts by means of integrals, thereby avoiding the need to match the wave function 
to known asymptotic functions at large distances. For two-body scattering calculations a 
representation of the T— matrix in momentum space is most convenient because this space 
is most closely related to the ingoing and outgoing momenta of the scattering process. 
However, for other applications it is convenient to formulate the T— matrix, and its respective 
integral equation, in configuration space. One such application occurs for the solution in 
configuration space of the three-body Faddeev integral equations, which require as input the 
two-body scattering matrices for each of the three arrangements [1]. 

The function T{E;r,r') satisfies the integral equation 

POO 

T{E; r, r') = V{r) 5{r - r') + V{r) / GoiE; r, f) T{E; f, r') df (1) 

Jo 

Here Qq is the undistorted (free) two-body Green's function Qq = {E + ie — Hq)~^, that, for 

a particular partial wave in configuration space with angular momentum i = is given by 

the well known expression 

GoiE- r, f) = —FiE- r<) G{E; r>) (2) 

where F{E;r) = sin{kr);G{E;r) = cos{kr) or exp(iA;r). In the case of £ ^ 0, F and G 
are replaced by the corresponding Bessel-Ricatti functions. Here k is the wave number in 
units of inverse length, and E and V{r) are the energy and the local two-body potential, in 
units of inverse length squared, respectively. In these units k = [Ey/'^ (for negative energies 

K={-Ef'^). 

The purpose of this paper is to present an algorithm for the solution of Eq. (1) for 
the two-variable function T{E;r,r') based on the use of spectral expansions in terms of 
Chebyshev polynomials [2], [3]. This method, denoted as I EM, has been found to be 
well suited for the solution of the one-variable Lippmann-Schwinger equation for the wave 
function ■?/'(r),both for positive [4], [5] and negative [6] energies. One essential feature of the 
lEM is to divide the entire radial range into partitions, and expand the unknown function 
in each partition in terms of a small number (~ 17) of basis functions such as Chebyshev 
or Legendre polynomials. Partitioning has the merit a) that the number of mesh points is 



very economical, since the partitions are concentrated automatically in those regions where 
the function varies most rapidly, b) the accuracy control of the spectral expansions is also 
maintained, and c) the singularity in the driving term, described below, can be easily handled 
by placing the boundary between two partitions at the position of the singularity. The basic 
reason that partitioning can still be maintained for the calculation of T is that the integral in 
Eq. (1) is over one variable only, i.e., the second variable is introduced only parametrically 
through the driving term V{r) 6{r — r'). 

The properties of T in configuration space are little known, mainly because T has usually 
been studied in the momentum representation. However, the potential between atoms are 
generally given in configuration space, and the presence of large repulsive cores as well as long 
range parts of the potential, that become important at ultra-cold energies for atom-atom 
scattering, can cause difficulties in the momentum representation. It is hoped that a good 
understanding of T in configuration space can provide a simpler way for overcoming these 
difficulties. Another important motivation is that the solution of the three-body Faddeev 
integral equations in momentum space [7], [8] achieves an accuracy of three to four significant 
figures [9], which is adequate for nuclear physics applications, but is inadequate for atomic 
physics applications. In order to make use of the lEM^s high accuracy for the solution of 
integral equations in configuration space, the three-body Faddeev integral equations have 
been formulated in configuration space [1]. Since a basic input into these equations are 
integrals over the two-body T— matrices [10], knowledge of the properties of the two-body 
T— matrix in configuration space provides an important contribution. 

II. THE FORMALISM 

An important relation between the T— matrix, the potential V, and the scattering wave 
function ip is 



V{r) ipir) = T{r,f) F{f)df. (3) 

Jo 

Here ip is the solution of the Schrodinger equation 

-^ + iV-k')^ = (4) 

or of the corresponding Lippmann- Schwinger (L-S) equation 

/•oo 

^(r) = F(r) + / goir,r') V{r') ^(r') dr'. (5) 

Jo 



Symbolically the solution of the L-S equation for ip can be written as ip = [1 — Qo V)~^F, 
from which one finds T = V {1 — Qo V)^^. The latter can be rewritten as T = (1 — VQo)~^V, 
from which either Eq. (1) follows, or the alternate form 

T{E;r,r') = V{r) 5{r-r') + V{r') T{E;r,f)go{E;f,r')df, (6) 

Jo 

obtained by using the formal identity (1 — VQo)^^y = V {1 — ^o^)^^- 

The function T{r, r') is symmetric in r, and r' as can be seen by comparing Eqs. (1) and 
(6) and making use of the symmetry of the Green's function Qq. Even though the T— matrix 
has a discontinuous derivative at r = r', as discussed below, the integral of the product of 
T with another function <l> 

f{E;r)= T{E,r,r') ^{r') dr' (7) 

Jo 

over one of the two variables is a continuous function of the other variable, and obeys a 
continuous one-variable integral equation 

POO 

f{E; r) = V{r) <l>(r) + V{r) / GoiE; r, f) T{E; f) df, (8) 

Jo 

that can be solved without difficulty [10] using methods described earlier [2]. 

The first step in solving Eqs. (1) or (6) is to separate out the delta function term from 
T by defining R{r, r') according to 

T{E; r, r') = V{r) 6{r - r') + R{E; r, r'). (9) 

As a result of Eq. (1) the L-S equation for R becomes 

R{E; r, r') = V{r) / GoiE; r, f) R{E; f, r') df + V{r)go{E; r, r')K(r'), (10) 

Jo 

which is the object of the present numerical evaluation. One physical interpretation of 
R{E; r, r') is that its integral over F represents all second and higher order Born terms of 
the iterative solution of the L-S Eq. (5) for the wave function 

^(r) = F{r) + / go{E;r,f) V{f)F{f) df 



JO 
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+ / go{E;r,f) / R{E;f,r') F{r')dr' df. (11) 

Jo Uo J 

This equation follows by inserting (1) into the right hand side of Eq. (3) and cancelling the 
potential V{r). The second term involving R in Eq. (11) is likely to cancel the first two 
terms in the repulsive core region of V since ip is small in this region, hence integrals over 
R alone are generally to be avoided. 



III. THE SPECTRAL METHOD 

An extension of the conventional spectral expansion method lEM [2] for solving Eq. (10) 
is described in this section. For simplicity the energy variable is suppressed in what follows, 
and the driving term V{r)QQ{E\ r, r')V{r') is denoted as D{E\ r, r'). In order to simplify the 
notation, it is useful to define the functions ^ and <& as 

^{E-r) = F{E-r)V{r) (12) 

(5{E-r) = G{E-r)V{r) (13) 



in terms of which 



D{E- r, r') = -\'S{E- r<) &{E- r>). (14) 



Since D has a discontinuous derivative at r = r', the solution of Eq. (10) also acquires the 
same type of discontinuity. 

Initially, without regard of the position of r', the radial interval is truncated at an upper 
limit < r < rmax , and that region is divided into M partitions. The size of each partition 
i, and the corresponding number M is determined automatically by the accuracy criterion 
according to the properties of Chebyshev expansions. [3], ([4]), and for this reason these 
partitions are denoted as Chebyshev partitions. The lower and upper limits of each partition 
are denoted by tj_i and U, respectively, 

U-i <ri< ti, (15) 

and if r is contained in partition i it is denoted as rj. A second mesh of Nrp equidistant 
points r' is set up in the radial interval [0, rmax] at a mesh distance of hrp. For a point 
r' on this mesh that falls into a particular Chebyshev partition j, this partition is further 
divided into two, with r' located at the border between the two sub-partitions. The left 
(right) sub-partition is denoted as Jl Ur)- 

In what follows the energy is no longer written explicitly, and the angular momentum is 
assumed to be zero. The integral operator Qo restricted to partition i is denoted as Qi and 
is defined as 

g,^ = -Icin) I ^ F{f)^{f) df - ^F{n) f " Gif)^if) df, (16) 



where ^ is an arbitrary function being operated upon by Qi. Inserting Eq. (9) into Eq.(lO) 
one obtains 

R^J - V{ri)G, Ri,, = D{r„r'j) + 0(r,) f ^ \-)^F{f) R{f,r'j)df 




rrmsLX 1 

+ ^{r^) {-)-G{f) R{f,r'^)df. (17) 

J ti 

where the left hand side of Eq. (17) denotes Rir'i^r^-) — Vir^i) f/ Qo(r^f) Rif^r') dr. The 
integrals in Eq.(17) do not depend on r, and are denoted as 

1 /-rmax 



Mrj) = -j^l G{f) R{f,r'^)df (18) 



1 



tj-i 



B,{r]) = —r F{f) R{f,r'j)df. (19) 

In view of the linearity of the operator (1 — VQi) on the left hand side of Eq. 17, the solution 
of Eq. 17 can be written as [2] a linear combination of three functions Y, Z, and S 

R (r„ r^) = A{r'^) F,(r,) + B,{r'^) Z,{r,) + S,,, (r„ r'^ (20) 

that are the solutions of 

[1 - V{n)g,] Y, = din) (21) 

[1 - V{r,)g,] Z, = (3{r,) (22) 

[1-V{n)g,] S^, = D{n,r'^) (23) 

in each interval i. Up to here the procedure is identical to the I EM for the solution of 
two-body L-S equations. A new feature is the appearance of the third function Sij, that, in 
view of Eq.(14) is given by 

S^J{n,r'^) = -^Y,in)i3ir'^) t<j (24) 

^.,(r„r;.) = -^^.(r,)^(r;) t>j. (25) 

For the partition j that subsequently was sub-divided into the left and right sub-partitions 
Jl and j/j, with r' at the border between the two partitions, Eqs. (24) and (25) remain valid 
for each of the sub-partitions, provided that i or j are replaced by the appropriate value Jl 
or JR. 



The calculations of the coefficients A and B is as follows. Inserting Eq.(20) into Eqs.(18) 
and (19) one obtains for i = 1,2, ..M 

M 

A^{r',) = Y. vn{r'.) (26) 

h=i+l 
i-l 

B^{r'^) = Y. Mr',), (27) 



h=l 



where 



Vhir'^) = -Ahir'.) (GF), - B^^ir'.) {GZ)r. - {GS)^,, (28) 

Uhir'j) = -A(rp (FF), - 5,(r;.) FZ)j, - {FS)h,j. (29) 



In the above, 



th ^ 



i^Oh = / TV{r)Ur) df, (30) 

where rj is either G or F , and ^h is either Yh or Zh, and 



^^^^h,j= TVir)Sh,jir,r'j) dr (31) 

From Eqs. (26) and (27) it follows that Bi = Am = 0. In view of Eqs. (24) and (25) one 
obtains for the inhomogeneous overlap terms 

{VS),, = -1{VY), (&{x,) h<j (32) 

and 

{vS)f,, = -I ivZ), dix,) h > J. (33) 

Here h stands for any of the partitions i = 1,2, ..M, with the exception than when h = j, h 
is replaced by either jl or j^. The equations above still apply, since r' is located either on 
the right or left border of the partition, respectively, as is explained in more detail below. 

A recursion relation between the coefficients A^ and B^ and those from a neighboring 
partition /i ± 1 can be obtained [4] by writing Eqs. Eqs. (26) and (27) in the form 

Ah+i -Ah = -Vh+i (34) 

Bh+i - Bh = Uh (35) 



According to Eqs. (28) and (29) the functions Vh and Uh contain the coefficients Ah and Bh-, 
and hence, after some algebra one obtains, for i + 1 and i both different from j [4] 

fi.+i I '_' I =rH ': I + I '^^'^^'•^ |;^ = 1,2,M-1, (36) 




i+l 




where 





^.+1 = (37) 



and 

'-(^^■>-^-<^^>-|- (38) 

1 

The matrices Vt and F have aheady been defined in Ref. [4], while the terms in Eq. (36) 
involving the functions S are new 

The method for obtaining the coefficients Ai and Bi as a function of r' consists in prop- 
agating the coefficients forward from partition i = 1 to partition ji, and backward, from 
partition M to ji- By equating the two results in partition j^, all the coefficients can be de- 
termined as a function of r'. These forward and backward propagations to the point r' bear 
some similarity to the method described in Ref. [6] for the calculation of energy eigenvalues 
of the Schrodinger equation. Further details for the calculation of the coefficients Ai and Bi 
are given in the Appendix 1. For i ^ j the result is a semi-separable expression of the form 

R (r„ r;.) = [A,(r;.) - &{r'^)/k] F,(r,) + B,(r;.) Z,(r,), i < j (39) 

R {n, r;.) = A{r'^) Y,{n) + [5,(r;) - d{r-)/k] Z,{n), t > j. (40) 

Integrals of R{r,r') over a known function $(r), as defined in Eq. (7), can be obtained 
by ffist calculating R and then performing the required integral, as is done below in order 
to check the accuracy of the results for R. However, a simpler and more accurate way is to 
solve the L-S equations (8) directly for these integrals. The present code can be modified 
for this purpose by replacing the driving term D in Eq. (14) by the continuous function 
V{r) $(r) (no right or left sub-partitioning being required), calculating the function S in each 
partition according to Eq. (23), choosing an appropriate point r' at the boundary between 
two Chebyshev partitions, performing the forward and back ward recursion relations for the 
quantities P,p and Q, q, respectively, and finally solving Eqs. (58) for the coefficients Ai and 



Bm- This operation has to be performed once only for the one point r' and can be carried 
out by the I EM with the same precision as the solution for the wave function [2], [4] or the 
binding energy [6]. 

If the driving term in Eq. (14) were set to zero, the quantities Sij in Eq. 23) would 
vanish, and hence Eq. (10) would be of the form 

R{E; r, r') = V{r) / GoiE; r, f) R{E; f, r') dr. (41) 

Jo 

The solution R of this equation vanishes, unless the determinant of the discretized form of 

the operator 1 — V{r) j^ df Qq{E\ r, f) is zero. This determinant does vanish for particular 

discrete energies which are the bound states of the two-body system. As the energy E ranges 

over the negative values which include bound state energies, poles in the i?-matrix occur. 

However, the exploration of these singularities will be left to a future study. 

1. Units and Dimensions 

The potential energy V and the energy E in the Schrodinger equation are normally given 
in units of energy, while the distance is in units of length {£). By multiplying all terms in the 
Schrodinger equation by 2/i//i^, where n is the reduced mass of the two colliding particles, 
and h is Planck's constant divided by 27r, then all terms in the scaled Schrodinger equation 
acquire dimensions of £~^.The new potential energy V and the energy E then are 

V = QV (42) 

E = k'^ = QE (43) 

where Q is a scaling factor. 

For nuclear physics applications E and V are given in MeV, the radial distance in fm 
and the corresponding value of Q is 

Qn = 2^1 /h\ (44) 

The corresponding scaled Schrodinger equation is Eq. (4), and the corresponding L-S equa- 
tion is (5). The dimension of the Green's function Qq., Eq. (2), is l^ (for the nuclear case 
I = fm) and hence the operator QQ{r,r')V{r')dr' has no dimension. In view of Eq. (3), 
and since F and ip have no dimension, the T— Matrix has dimension i~^. This dimension 

9 



is compatible with Eq. (9), since the delta function has dimension i ^ . The dimension of R 
is also £~^, which is compatible with Eq. (10). The dimension of A and B, Eqs. (18) and 
(19), is £-\ the dimension of Y and Z, Eqs. (21)-(22), is r^, and that of S, (23) is r^. 
The quantities {r]^) , Eq. (30), and {i]S) Eq. (31), have dimensions i^ and £^^, respectively. 
The quantities or ^ have dimensions i^^, and hence Eqs. (39) and (40) are dimensionally 
self-consistent, each term having the dimension i~^. 

For atomic physics applications E is given in atomic energy units 2R, (M ~13.606 eV) 
and the distance r in units of the Bohr radius ao, in which case 

Qa = %l X 2M = ^, (45) 

h^ me 

where rUe is the mass of the electron. In this case E and r, and similarly the scaled quantities 
V and E, Eqs. (42) and (43), respectively, have no dimension and hence all the quantities 
described above for the nuclear case have no dimensions either. 

IV. NUMERICAL EXAMPLES 

Two different potentials in Eq. (4) are used for the numerical examples presented here. 
One is a simple exponential potential and the other is a potential describing the interaction 
between two Helium atoms, based on the potential TTY [11], [12]. For the He-He case a soft 
repulsive core is applied for a distance less than Vcut = 4.5. The numerical lEM calculations 
are carried out from r = to r = r^-a^x, so as to include the effect of the core accurately. 
The value of the accuracy parameter is tol = 10~^, and the wave number is k = 1.5. The 
potential and the wave number are in units of inverse length squared, and inverse length, 
respectively. 

A. The exponential example. 

The potential in Eq. (4) is V{r) = ±exp(— r). The i?-matrix for the repulsive case is 
shown in Fig (1), and the behavior of the diagonal part of R is shown in Fig. (2). The 
derivative discontinuity of R{r,r') at r = r' is clearly visible in Fig. (1), and it is also seen 
that R is large where the potential is large. These features are also present for the He — He 
case, described below. 
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r prime 







FIG. 1: The i?-matrix for the exponential potential V = exp(— r). The wave number k isl.5. The 
discontinuity in the derivative, for r = r' is clearly visible. 
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FIG. 2: The diagonal part of the i?— matrix, shown in Fig. (1). 
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FIG. 3: The He — He atom interaction, based on the potential TTY [11]. The units for the 
distance are ag (TheBohr radius) and (ao)^^ for the potential. For distances less than 4.5 gq a 
repulsive soft core is matched to this potential, as described in the text. At r ~ its value is ~ 150 
(ao)~^ . The He — He dimer is very weakly bound, with a binding energy of 2.58 x 10~^ (oq)"^, 
or ~ 10^'^ey. Upon division by the factor 7296.3 this potential is transformed to atomic energy 
units, as described in Ref.. [6]. 

B. The He-He case. 

For the He — He case the potential is shown in Fig. (3) 

It has a strong repulsive core, that is not adequately represented by the formula described 
in Refs. [11], [12]. Since the result for the i?— matrix depends significantly on the nature of 
the repulsive core, an "artificial" repulsive soft core was used to replace the repulsive core 
of TTY for distances < r < rcut- This core is of the form a + br + cr'^, with the coefficients 
a, 6,and c determined such that the repulsive core is equal to the value, the first, and second 
derivatives of TTY at r = rcut- This form of the core is less unphysical than the "cut-ofF' 
core used in our previous numerical calculations of the binding energy of the He — He dimer 
[6], for which V was replaced by a constant value V{rcut) for < r < Vcut- 



12 




I -400 



r prime 



0.4 0.4 



FIG. 4: The i?— matrix for the He — He potential at distances weh inside of the repulsive core, 
with rcut = 4.5, and k = 1.5. 

The i?— matrix at small distances is shown in Fig. (4). The values of R are very large in 
the region or the large repulsive core, but they become smaller at larger distances. Figs. (5) 
and (6), for which the potential becomes progressively smaller. 

The diagonal parts of the He — He /?— matrix are shown in Figs. (7) and (8), for short 
and large distances, respectively. 

C. Accuracy Tests. 



The accuracy of the functions Y and Z in each partition, Eqs. (21) and (22), determined 
by the size of the partitions, is better than tol = 10~^. In order to obtain a measure of the 
loss of accuracy that takes place in the subsequent steps of the calculation, the following 
two tests are performed: one that checks the symmetry of the i?— matrix, and the other that 
performs the numerical integral 



3(t) 



R{r,r') X F{r)dr 
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FIG. 5: The He-He i?— matrix, same as in Fig. (4), but at larger distances, straddling the end of 
the repulsive core and the attractive potential valley. 




r prime 



9 9 



FIG. 6: The He- He i?— matrix, same as in Fig. (4), at distances between 9 and 12 a^ where the 
potential is of the r~^ form. The value of R is of the order of 10"^ {ao)~^. 
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FIG. 7: The diagonal part of the He — He i?-matrix, shown in Fig. (4). at sniah distances 
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FIG. 8: Same as Fig. (7) for large distances. 
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n 


^ _ nXTT 
' 16 




Exp. potil 


He - He pot'l 


10 


-1.96 


L 


3.85431605 [ - 2] 


-1.079593140 [1] 






R 


3.85431605 [-2] 


-1.079593140 [1] 


15 


~2.95 


L 


4.05244938 [-3] 


2.34862829 [1] 






R 


4.05244931 [-3] 


2.34862833 [1] 


20 


-3.93 


L 


-5.633010764 [-3] 


4.12156423 






R 


-5.633010723 [-3] 


4.12156414 


25 


~4.91 


L 


-1.073578042 [-3] 


1.4144533 [ - 2] 






R 


-1.073578029 [-3] 


1.4144530 [-2] 


30 


-5.89 


L 


7.19293011 [-4] 


8.2459426 [ - 2] 






R 


7.19293028 [-4] 


8.2459429 [-2] 


35 


-6.87 


L 


2.05522067 [-4] 


3.80032458 [-2] 






R 


2.05522071 [-4] 


3.80032395 [-2] 


40 


-7.85 


L 


-8.59613654 [-5] 


1.9356054 [ - 2] 






R 


-8.59613675 [-5] 


1.9356055 [-2] 



TABLE I: Values of the right (R) and left (L) sides of Eq. (46). Numbers in square brackets denote 
powers of 10 



which, according to Eq. (3), should equal [ip{r') — F{r')] V{r') 

J(r')=? lijir') - F{r')] \/(r'), 



(46) 



In Table I the values of the two sides of the above equation are compared in order to check 
the accuracy of their agreement with each other. For the exponential case the number of 
Chebyshev support points in each partition is A^ = 17, and only four partitions are required 
to cover the interval from r = to r = 25. 

For He — He case, if iV = 17, the interval from to 250 requires 26 partitions, 11 of 
which lie in the region of the repulsive core. The corresponding accuracy of the integral 3{r) 
(not displayed in the table) is three to four significant figures. As A^ is increased, the size of 
each partition increases, correspondingly the number of partitions decreases. The accuracy 
of R increases the smaller the number of partitions, because the Chebyshev support points 
are the ones that carry the large change in the i?— values, rather than the coefficients A and 
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B.. The results for 3{r) for A^ = 65 are shown in "L" hnes of Table (I), from which one can 
deduce an accuracy for R of at least 7 significant figures. For this case there are only four 
partitions in the region [0 — > 5] of the repulsive core, and four more partitions that cover 
the whole remaining distance [5 — > 250]. These accuracy results for both the exponential 
and the He — He cases are confirmed by the degree of symmetry of the i?— matrix. 

Using MATLAB, version 7.0.1.24704(i?14) on a 2.8 GHz Intel computer, the total time 
required for a He — He calculation of R{r,r'), with n = 1273 equidistant meshpoints for 
each r and r', rmax = 250, A^ = 17, and tol = 10^^ requires 114 seconds. (This time scales 
like n^) For the same calculation, but with A^ = 65, the total MATLAB time is 261 s. The 
corresponding calculation of the wave function ip, for r^ax = 1500 and also 1273 equidistant 
meshpoints takes between 1 and 2 seconds. 

V. SUMMARY AND CONCLUSIONS 

In this paper the spectral lEM method for solving the one variable Lippmann-Schwinger 
integral equation {L — S) in configuration space for the wave function [2], [4] is extended to 
the calculation of the two variable two-body T— matrix, also in configuration space. The 
T-matrix is written as the sum of a delta-function multiplying the (local) potential plus 
a reminder, called the i?-matrix. The L — S equation for the latter is the object of the 
calculation. The main complication is that the driving term for this equation is continuous 
but has a discontinuous derivative at r = r'. That discontinuity also propagates into the 
i?-matrix, and makes the calculation more cumbersome. This difficulty is overcome by 
constructing the partitions into which the radial domain is divided in such a way that the 
point r' is located at the end-point of a partition. This procedure is repeated for all values 
of r' contained in an equidistant set of mesh of points. Otherwise the calculation is very 
similar to the spectral I EM procedure for the solution of the one variable L — S equation, 
whose accuracy and reliability was previously demonstrates in various applications [5], [6]. 
The present method complements a previous investigation of expanding T in terms of a sum 
over Sturmian functions [13]. 

Two numerical examples are described. One, in which the potential is a simple exponen- 
tial function of r, and the other in which the potential describes the interaction between two 
Helium atoms [11]. For the exponential case the high accuracy, which is determined by the 
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inputted tolerance parameter tol, is maintained. For tol = 10~*, the i?— matrix is accurate 
to 7 significant figures. For the He — He case the accuracy of R is reduced to between 3 and 
4 significant figures. This loss of accuracy is attributed to the presence of the large repulsive 
core present in the He — He interaction. If the number of Chebyshev support points in each 
partition is increased from 17 to 65, then the accuracy of the i?— matrix for the He — He 
case increases to 7 significant figures. 

The calculation of integrals T{r) over T{r, r') can be carried [10] out with far less com- 
putational effort than the calculation of T{r.r'), since these integrals can be obtained from 
the solution of a L — S* equation (8) that is very similar to the L — S equation for the wave 
function, whose accuracy was found to be close to 7 or 8 significant figures [6]. This is 
of interest, since the integral Faddeev equations for the three-body system, formulated in 
configuration space [1] have such integrals as important inputs. It is our aim to investigate 
the numerical solution of the three-body integral Faddeev equations in configuration space 
in a future study, using the results of the present paper, and that of Ref. [10]. 

Acknowledgments: The author is greatly indebted to Professor Walter Glockle 
for stimulating and consistent help with this work, and for pointing out its underlying 
importance with regards to the envisaged three-body calculations. The author is also grateful 
to Professor I. Koltracht, at UConn, for help with the computational-mathemathical aspects 
of this work. 

Appendix 1. The calculation of the coefficients A and 
B. 

Implementing the solution of the recursion relations (36) in a numerically efficient way 
will be described below. A short summary is as follows. By implementing the forward 
recursion relations (47) described below, starting with {Ai,Bi) = (1,0) one obtains the 
quantities Pi and pi, Eq. 48), in terms of which the corresponding values of Ai and Bi for 
i < j, are given by {Ai,Bi)'^ = AiPi+ pi, Eq. (56) with Ai still unknown (here T means 
"transposed"). Similarly, the coefficients Ai and Bi are propagated backwards from the 
last partition M also into ji, using Eq. (51), starting with {Am,Bm) = (0, 1) , Eq. (56). 
The result is given by a combination of the quantities Qi and g^, Eq. 52). In terms of 
these quantities the corresponding coefficients Ai and Bi for and for i > j are given by 
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{Ai,Bi)'^ = BMQi+ Qi, Eq. (57) with Bm still unknown. By equating the coefficients Aj^^ 
and Bjj^ in the same partition jl, obtained by the two procedures, one obtains two equations 
for the coefficients Ai and Bm in terms of the known "driving" terms Pi and g^ involving 
the functions Sij. These equations (58) then can be solved for Ai and Bm, but since they 
depend on the location of r', the resulting values of Ai and Bm also depend on r'. By 
repeating this procedure for all values of r' in partition j for all partitions j one obtains Ai 
and Bi for all values of r' . 

In detail the procedure is as follows. In view of Eq. (36) the "forward" recursion relations 
are 




{n 



i + ly 



r. 




+ 



e(r': 



k 




1,2,.. J -2. (47) 



In the above, when i = j — 1 then i + 1 = Jl- By repeatedly using Eqs. (47) one can express 
Ai and Bi in terms of Ai and Bi, with the result 



A-f 
Bf 



Pi+PiiJ), i = 2,S,..J -1,JL, 



(48) 



where the recursion for the Pi and pi column vectors for i = 1,2, ..j — 1 are 



Pi 



-* j+l — * 'i+l J- j -* 15 



(49) 



Pi 







; Pi+iir'j) = ^i+i 



T,pM) + -eirr 



k ' ^ 




When i = j — 1, the above expressions yield Pj^ and pj^ 
The "backward" recurrence relations are 




(r 



\-i 



j-ly 



n. 



A" 
B' 



j-i 






k 




(50) 



i = M,M-l,..,j+l. (51) 



When i = j + 1 then i — 1 is replaced by jji. The superscripts / and b stand for "forward" 
and "backward", respectively. By repeatedly using Eq. (51) one obtains 



A' 
B' 



Qi + qr{j), i = M-i,..,j + i,jR 



(52) 
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and the corresponding recursion relations for the Qi and qi column vectors are 



Ql — \ I 5 Qi-l — Tj-l ^i Qi 




Qi 




; Qi^iirj, 



j-i 



^i Qi{r'j) + 



k 




(53) 



(54) 



In these expressions r' stands to the left of any partition, hence Eqs. (51) apply, and 
Eqs. (53) are valid for i — 1 = M — 2, M — 1, .., Jl while Eqs. (54) are valid only for 
i — 1 = M — 2,M — 1, ..,JR. For the backward recursion from jr to Jl the expression for g^^ 
has to be replaced by 



<l,Ar')=T 



]L 



^hn. IJRVj 




(55) 



^(-^•) (GZh, 

The vectors Pi for i < j and Qi for i > j are independent of r' and the corresponding 
vectors Pi and Qi all dependent on r' via the factor (3{rj)/k and '^{r'j)/k respectively, a 
feature that provides economy in the numerical calculation. The vectors {A^ ^ B^)'^ as given 
by Eq. (48) are the values of the coefficients A and B when the values of (Ai, Bi)^ in the 
first partition are set to (1, 0)^ and similarly for (A^, -B^)"^, as given by Eq. (52), when the 
values of {Am-, BmY^ in the last partition M are set to (0, 1)"^. The correct values of (A, B)^ 
are obtained by multiplying the vectors Pi by Ai and the vectors Qi by -Bm, i-e., 

Ai(r;.)P.+p.(r;.), ^ = 2,3,..,j-l,ji, (56) 




BM{r[.) Qi + g,(r'), ^ = M - 1, .., j + 1, j^. 



(57) 



In order to determine the values of Ai and Bm one equates to each other the values of 
{A^B)^ obtained by Eqs. (56) and (57), respectively, in the partition j^ with the result 



p(l) _Q{1) 
p(2) _Q(2) 



]L 




QjL - PjL 



(58) 



The procedure described above is repeated for every value of r' in partition j. A convenient 
set of such values, denoted as r' ^i, s' = 1, 2, ..A^ + 1, which originate from in the equidistant 
mesh of points r' and are located in the Chebyshev mesh j. 



20 



The case when i = j is as follows, but for clarity of notation the subscript j, is suppressed. 
Each Ls' and Rs' partition has its own A^ + 1 Chebyshev support points by means of which 
the Chebyshev expansion coefficients of the functions Y^ , and Z^ ,, as well as Yr , and 
Z^ , are calculated via the I EM in the Left and Right sub-partitions. By means of these 
expansion coefficients the values of the functions Y and Z, and hence also the values of these 
functions at a point r contained in the equidistant set of mesh-points can be calculated, and 
hence the functions R can be obtained at the equidistant mesh points of r and r' . In a 
simplified notation, Rj{s,s') = [A{s') — (3{s')/k] Yl_^,{s) + B{s') Zl^,{s), s < s', and a 
similar expression based on Eq. (40) for s > s'. 

After the cases that the points r are contained in different Chebyshev partitions than 
point r', are combined with the case where r and r' are contained in the same Chebyshev 
partition, one obtains the final result for R{rs,r'^,). The dimension of this matrix is equal to 
the number Nrp of equidistant mesh points. 
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